!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \brief Definition of mathematical constants and functions.
!> \par History
!>      Adapted for use in CP2K (JGH)
!>      JGH (16-06-2002) : Added Gamma functions
!>      JGH (10-08-2004) : Added Euler constant (gamma)
!> \author Matthias Krack
! **************************************************************************************************
MODULE mathconstants

   USE kinds,                           ONLY: dp

   IMPLICIT NONE

   PRIVATE

   PUBLIC :: maxfac, pi, twopi, fourpi, rootpi, oorootpi, zero, one, &
             degree, radians, gaussi, fac, ifac, dfac, gamma0, gamma1, euler, &
             sqrthalf, sqrt2, sqrt3, sqrt5, sqrt7, sqrt15, sqrt21, sqrt35, &
             sqrt105, z_mone, z_one, z_zero

   ! Factorial function fac
   ! Inverse Factorial function ifac
   ! Double factorial function dfac
   ! Gamma functions
   ! gamma(n) = gamma0(n) = (n - 1)!
   ! gamma(n + 1/2) = gamma1(n) = SQRT(pi)/2^n (2n - 1)!!

   INTEGER, PARAMETER :: maxfac = 30
   REAL(KIND=dp), PARAMETER, DIMENSION(0:maxfac) :: fac = (/ &
                                      0.10000000000000000000E+01_dp, 0.10000000000000000000E+01_dp, 0.20000000000000000000E+01_dp, &
                                      0.60000000000000000000E+01_dp, 0.24000000000000000000E+02_dp, 0.12000000000000000000E+03_dp, &
                                      0.72000000000000000000E+03_dp, 0.50400000000000000000E+04_dp, 0.40320000000000000000E+05_dp, &
                                      0.36288000000000000000E+06_dp, 0.36288000000000000000E+07_dp, 0.39916800000000000000E+08_dp, &
                                      0.47900160000000000000E+09_dp, 0.62270208000000000000E+10_dp, 0.87178291200000000000E+11_dp, &
                                      0.13076743680000000000E+13_dp, 0.20922789888000000000E+14_dp, 0.35568742809600000000E+15_dp, &
                                      0.64023737057280000000E+16_dp, 0.12164510040883200000E+18_dp, 0.24329020081766400000E+19_dp, &
                                      0.51090942171709440000E+20_dp, 0.11240007277776076800E+22_dp, 0.25852016738884976640E+23_dp, &
                                      0.62044840173323943936E+24_dp, 0.15511210043330985984E+26_dp, 0.40329146112660563558E+27_dp, &
                                      0.10888869450418352161E+29_dp, 0.30488834461171386050E+30_dp, 0.88417619937397019545E+31_dp, &
                                                    0.26525285981219105864E+33_dp/)
   REAL(KIND=dp), PARAMETER, DIMENSION(0:maxfac) :: ifac = (/ &
                                      0.10000000000000000000E+01_dp, 0.10000000000000000000E+01_dp, 0.50000000000000000000E+00_dp, &
                                      0.16666666666666666667E+00_dp, 0.41666666666666666667E-01_dp, 0.83333333333333333333E-02_dp, &
                                      0.13888888888888888889E-02_dp, 0.19841269841269841270E-03_dp, 0.24801587301587301587E-04_dp, &
                                      0.27557319223985890653E-05_dp, 0.27557319223985890653E-06_dp, 0.25052108385441718775E-07_dp, &
                                      0.20876756987868098979E-08_dp, 0.16059043836821614599E-09_dp, 0.11470745597729724714E-10_dp, &
                                      0.76471637318198164759E-12_dp, 0.47794773323873852974E-13_dp, 0.28114572543455207632E-14_dp, &
                                      0.15619206968586226462E-15_dp, 0.82206352466243297170E-17_dp, 0.41103176233121648585E-18_dp, &
                                      0.19572941063391261231E-19_dp, 0.88967913924505732867E-21_dp, 0.38681701706306840377E-22_dp, &
                                      0.16117375710961183490E-23_dp, 0.64469502843844733962E-25_dp, 0.24795962632247974601E-26_dp, &
                                      0.91836898637955461484E-28_dp, 0.32798892370698379102E-29_dp, 0.11309962886447716932E-30_dp, &
                                                    0.37699876288159056439E-32_dp/)
   REAL(KIND=dp), PARAMETER, DIMENSION(-1:2*maxfac + 1) :: dfac = (/ &
                                      0.10000000000000000000E+01_dp, 0.10000000000000000000E+01_dp, 0.10000000000000000000E+01_dp, &
                                      0.20000000000000000000E+01_dp, 0.30000000000000000000E+01_dp, 0.80000000000000000000E+01_dp, &
                                      0.15000000000000000000E+02_dp, 0.48000000000000000000E+02_dp, 0.10500000000000000000E+03_dp, &
                                      0.38400000000000000000E+03_dp, 0.94500000000000000000E+03_dp, 0.38400000000000000000E+04_dp, &
                                      0.10395000000000000000E+05_dp, 0.46080000000000000000E+05_dp, 0.13513500000000000000E+06_dp, &
                                      0.64512000000000000000E+06_dp, 0.20270250000000000000E+07_dp, 0.10321920000000000000E+08_dp, &
                                      0.34459425000000000000E+08_dp, 0.18579456000000000000E+09_dp, 0.65472907500000000000E+09_dp, &
                                      0.37158912000000000000E+10_dp, 0.13749310575000000000E+11_dp, 0.81749606400000000000E+11_dp, &
                                      0.31623414322500000000E+12_dp, 0.19619905536000000000E+13_dp, 0.79058535806250000000E+13_dp, &
                                      0.51011754393600000000E+14_dp, 0.21345804667687500000E+15_dp, 0.14283291230208000000E+16_dp, &
                                      0.61902833536293750000E+16_dp, 0.42849873690624000000E+17_dp, 0.19189878396251062500E+18_dp, &
                                      0.13711959580999680000E+19_dp, 0.63326598707628506250E+19_dp, 0.46620662575398912000E+20_dp, &
                                      0.22164309547669977187E+21_dp, 0.16783438527143608320E+22_dp, 0.82007945326378915594E+22_dp, &
                                      0.63777066403145711616E+23_dp, 0.31983098677287777082E+24_dp, 0.25510826561258284646E+25_dp, &
                                      0.13113070457687988603E+26_dp, 0.10714547155728479551E+27_dp, 0.56386202968058350995E+27_dp, &
                                      0.47144007485205310027E+28_dp, 0.25373791335626257948E+29_dp, 0.21686243443194442612E+30_dp, &
                                      0.11925681927744341235E+31_dp, 0.10409396852733332454E+32_dp, 0.58435841445947272053E+32_dp, &
                                      0.52046984263666662269E+33_dp, 0.29802279137433108747E+34_dp, 0.27064431817106664380E+35_dp, &
                                      0.15795207942839547636E+36_dp, 0.14614793181237598765E+37_dp, 0.86873643685617511998E+37_dp, &
                                      0.81842841814930553085E+38_dp, 0.49517976900801981839E+39_dp, 0.47468848252659720789E+40_dp, &
                                       0.29215606371473169285E+41_dp, 0.28481308951595832474E+42_dp, 0.17821519886598633264E+43_dp/)
   REAL(KIND=dp), PARAMETER, DIMENSION(0:maxfac) :: gamma0 = (/ &
                                      0.00000000000000000000E+00_dp, 0.10000000000000000000E+01_dp, 0.10000000000000000000E+01_dp, &
                                      0.20000000000000000000E+01_dp, 0.60000000000000000000E+01_dp, 0.24000000000000000000E+02_dp, &
                                      0.12000000000000000000E+03_dp, 0.72000000000000000000E+03_dp, 0.50400000000000000000E+04_dp, &
                                      0.40320000000000000000E+05_dp, 0.36288000000000000000E+06_dp, 0.36288000000000000000E+07_dp, &
                                      0.39916800000000000000E+08_dp, 0.47900160000000000000E+09_dp, 0.62270208000000000000E+10_dp, &
                                      0.87178291200000000000E+11_dp, 0.13076743680000000000E+13_dp, 0.20922789888000000000E+14_dp, &
                                      0.35568742809600000000E+15_dp, 0.64023737057280000000E+16_dp, 0.12164510040883200000E+18_dp, &
                                      0.24329020081766400000E+19_dp, 0.51090942171709440000E+20_dp, 0.11240007277776076800E+22_dp, &
                                      0.25852016738884976640E+23_dp, 0.62044840173323943936E+24_dp, 0.15511210043330985984E+26_dp, &
                                      0.40329146112660563558E+27_dp, 0.10888869450418352161E+29_dp, 0.30488834461171386050E+30_dp, &
                                                    0.88417619937397019545E+31_dp/)
   REAL(KIND=dp), PARAMETER, DIMENSION(0:maxfac) :: gamma1 = (/ &
                                      0.17724538509055160273E+01_dp, 0.88622692545275801365E+00_dp, 0.13293403881791370205E+01_dp, &
                                      0.33233509704478425512E+01_dp, 0.11631728396567448929E+02_dp, 0.52342777784553520181E+02_dp, &
                                      0.28788527781504436100E+03_dp, 0.18712543057977883465E+04_dp, 0.14034407293483412599E+05_dp, &
                                      0.11929246199460900709E+06_dp, 0.11332783889487855673E+07_dp, 0.11899423083962248457E+08_dp, &
                                      0.13684336546556585726E+09_dp, 0.17105420683195732157E+10_dp, 0.23092317922314238412E+11_dp, &
                                      0.33483860987355645697E+12_dp, 0.51899984530401250831E+13_dp, 0.85634974475162063871E+14_dp, &
                                      0.14986120533153361177E+16_dp, 0.27724322986333718178E+17_dp, 0.54062429823350750447E+18_dp, &
                                      0.11082798113786903842E+20_dp, 0.23828015944641843260E+21_dp, 0.53613035875444147334E+22_dp, &
                                      0.12599063430729374624E+24_dp, 0.30867705405286967828E+25_dp, 0.78712648783481767961E+26_dp, &
                                      0.20858851927622668510E+28_dp, 0.57361842800962338401E+29_dp, 0.16348125198274266444E+31_dp, &
                                                    0.48226969334909086011E+32_dp/)

   ! Constants related to Pi

   REAL(KIND=dp), PARAMETER :: pi = 3.14159265358979323846264338_dp ! Pi
   REAL(KIND=dp), PARAMETER :: pio2 = 1.57079632679489661923132169_dp ! Pi/2
   REAL(KIND=dp), PARAMETER :: twopi = 6.28318530717958647692528677_dp ! 2*Pi
   REAL(KIND=dp), PARAMETER :: fourpi = 12.56637061435917295385057353_dp ! 4*Pi
   REAL(KIND=dp), PARAMETER :: rootpi = 1.77245385090551602729816748_dp ! SQRT(Pi)
   REAL(KIND=dp), PARAMETER :: oorootpi = 0.56418958354775628694807945_dp ! 1/SQRT(Pi)

   ! Euler's constant (Euler-Mascheroni constant)

   REAL(KIND=dp), PARAMETER :: euler = 0.57721566490153286060651209_dp

   ! Trivial constants

   REAL(KIND=dp), PARAMETER :: zero = 0.0_dp, &
                               half = 0.5_dp, &
                               one = 1.0_dp

   REAL(KIND=dp), PARAMETER :: sqrthalf = 0.70710678118654752440084436_dp, & ! SQRT(1/2)
                               sqrt2 = 1.41421356237309504880168872_dp, & ! SQRT(2)
                               sqrt3 = 1.73205080756887729352744634_dp, & ! SQRT(3)
                               sqrt5 = 2.23606797749978969640917367_dp, & ! SQRT(5)
                               sqrt7 = 2.64575131106459059050161575_dp, & ! SQRT(7)
                               sqrt15 = 3.87298334620741688517926540_dp, & ! SQRT(15)
                               sqrt21 = 4.58257569495584000658804719_dp, & ! SQRT(21)
                               sqrt35 = 5.91607978309961604256732829_dp, & ! SQRT(35)
                               sqrt105 = 10.24695076595959838322103868_dp ! SQRT(105)

   ! Conversion factors

   REAL(KIND=dp), PARAMETER :: degree = 180.0_dp/pi, & ! radians -> degree
                               radians = one/degree ! degree  -> radians

   COMPLEX(KIND=dp), PARAMETER :: gaussi = (0.0_dp, 1.0_dp) ! i = SQRT(-1)
   COMPLEX(KIND=dp), PARAMETER :: z_mone = (-1.0_dp, 0.0_dp), &
                                  z_one = (1.0_dp, 0.0_dp), &
                                  z_zero = (0.0_dp, 0.0_dp)

END MODULE mathconstants
